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ABSTRACT 

We propose a time-dependent one-zone model based on a quasi-Maxwellian 
'pile-up' distribution in order to explain the time-averaged high energy emission 
of TeV blazars. The instantaneous spectra are the result of the synchrotron 
and synchrotron self-Compton emission (SSC) of ultra-relativistic leptons. The 
particle energy distribution function (EDF) is computed in a self-consistent way, 
taking into account an injection term of fresh particles, a possible pair creation 
term, and the particles radiative cooling. The source term is not a usual power- 
law but rather a 'pile-up' distribution, which can result from the combination 
of a stochastic heating via second order Fermi process and radiative cooling. 
To validate this approach, we have performed time-averaged fits of the well- 
known TeV emitter Mrk 501 during the 1997 flaring activity period taking into 
account the attenuation of the high energy component by cosmic diffuse infrared 
background (DIrB) and intrinsic absorption via the pair creation process. The 
model can reproduce very satisfactorily the observed spectral energy distribution 
(SED). A high Lorentz factor is required to avoid strong pair production; in the 
case of smaller Lorentz factor, an intense flare in the GeV range is predicted due 
to the sudden increase of soft photons density below the Klein-Nishina threshold. 
The possible relevance of such a scenario is discussed. 
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1. Introduction 

It is now widely admitted that radio-loud active galactic nuclei (AGN) harbor mag- 
netized accretion-ejection structures involving a supermassive black-hole as a central en- 
gine. The EGRET experiment aboard the Compton-GRO satellite discovered more than 80 
gamma-ray emitting AGNs, all of them belonging of the blazar class (non-thermal contin- 
uum spectrum, optical polarization, flat radio spectrum and strong variability in all frequency 
bands). Some of these objects have been also firmly detected by Atmospheric Cerenkov Tele- 
scope (ACT) with an emission above 1 TeV. The two prototypes of TeV blazars are Mrk 
421 (Punch et al. 1992) and Mrk 501 (Quinn et al. 1996), two objects relatively close to 
us and roughly at the same distance, respectively z s ~ 0.031 and z s ~ 0.034. Thanks to 
the development of the ground-based gamma-ray astronomy the sample of the TeV emitters 
is increasing. During the last decade, several ACT teams have reported the detection or 
the confirmation of new sources : 1ES 1426+428 (Horan et al. 2002; Djannati-Ata'i et al. 
2002; Aharonian et al. 2002), 1ES 1959+650 (Nishiyama et al. 2000; Aharonian et al. 2003; 
Holder et al. 2003), 1ES 2344+514 (Catanese et al. 1998) and PKS 2155-304 (Chadwick et al. 
1999). A characteristic feature of blazars is the strong non-thermal emission from the radio 
to the gamma-ray range attributed to a relativistic jet supposed to be closely aligned with 
the observer light-of-sight. Their spectral energy distribution is quite typical and consists in 
two broad bumps. In the context of the Synchrotron Self-Compton (SSC) models, the low 
energy component, peaking in the X-ray domain for TeV blazars, is commonly attributed 
to the synchrotron emission of ultra-relativistic particles plunged into magnetic field. The 
second one is thought to be the result of up-scattering of the synchrotron photon field by 
the same population of ultra-relativistic particles via the inverse Compton (IC) mechanism 
(Jones, O'dell, & Stein 1974; Konigl 1981; Ghisellini, Maraschi, & Treves 1985). This model 
gives a good framework to explain the correlated variability for the the high and the low 
energy components. Even if the spectral properties of these objects seem to be understood, 
the different models do not discuss the origin and the physical mechanism of particle accel- 
eration. To reproduce the curved shape of the synchrotron and IC spectra on a wide energy 
domain, several authors have chosen a particle EDF parameterized by a simple or a broken 
power-law on a prescribed energy range [7min, 7max]- This choice is purely phenomenologi- 
cal and have no theoretical justification, even if in some special cases of shock acceleration 
(first order Fermi process) power- law EDF are expected (Jones 1994). For example, to 
reproduce the X-ray synchrotron bump, several authors use a simple power-law 71(7) oc 7~ s , 
7 £ [7min, 7max] but the dynamical range i.e. the ratio of 7max/7min is less than 10 (Pian et 
al. 1998). In this case, it seems to be more appropriate to consider a quasi mono-energetic 
distribution. In this work, we propose another primary type of EDF for emitting particles in 
order to reproduce the peculiar spectral energy distribution of TeV blazars. We assume that 
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the acceleration mechanism combined to radiative losses or/and an escape process produces 
a quasi-Maxwellian or 'pile-up' distribution, which is injected in a spherical region where it 
cools freely. The effect of cooling is to produce naturally a E~ 2 power-law in some limited 
range of energy. We also take into account the time dependence of the EDF to compare with 
the observations, considering that the observed spectra are always time-averaged spectra of 
intrinsically highly variable objects. In section 2, we present our kinetic scenario to obtain 
the energy spectrum of the particles in a self-consistent way, and we shortly describe the 
emission processes used to reproduce the blazar spectra. Finally we illustrate our approach 
in section 3, giving some results of SED fitting before concluding. 



2. The model 

2.1. Stochastic particles acceleration 

In the following, we will consider only a homogeneous one- zone model where all physical 
quantities are assumed to be averaged over the volume of the emission region. All spatial 
dependences are dropped from the equations. The particles distribution function f(p;t) is 
assumed to be isotropic in some frame, called the 'blob frame', moving relativistically with 
a bulk Lorentz factor r#. In this frame, it depends only on the modulus of the momentum 
p = \p\ and the time t. For relativistic particles the energy is given by E = r ym c c 2 ~ pc and 
the differential number density of pairs ^±(7; t) of reduced energy 7 is related to EDF 7(7) by 
the usual relation (time is implicit) dn± = n±{ r y) fry = 4iip 2 f(p) dp ~ 47r(m c c 2 ) 3 7 2 /(7) d7. 

We assume that the particles are accelerated stochastically by energy exchanges with 
resonant plasma waves in a weak turbulent medium. In our model, the acceleration zone must 
be localized : it could be the basis of a jet, or localized reconnection sites, or the interface 
between a relativistic beam and a confining jet as proposed for example by Henri & Pelletier 
(1991) in the framework of the 'two-flow model' (Pelletier 1985; Pelletier & Sol 1992). This 
insures that the particles will spend only a tiny fraction of time in the acceleration zone, be- 
fore being injected in a larger region where they cool freely. According to quasi-linear theory, 
the acceleration process can be described by a diffusion equation in the momentum space 
leading to a Fokker-Planck equation. This equation gives the time-dependent evolution of 
any initial particle density submitted to deterministic continuous energy changes or diffusive 
Markovian processes. We suppose that the characteristic acceleration time-scale is short 
compared with the other time-scales in the problem and we will focus our attention onto the 
stationary solution 7(7) of the Fokker-Planck equation. The diffusion coefficient ^ 77 (7) in 
phase space can be chosen as a power-law in terms of the Lorentz factor 7 (Lacombe 1977; 
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Henri & Pelletier 1991; Dermer, Miller, & Li 1996) 



where r e [1, 2] is the index of wave turbulent spectrum , assumed to be itself a power-law 
(e.g. r = 5/3 for a Kolmogorov turbulence, r = 3/2 for Kraishman one). The steady-state 
differential energy spectrum resulting from a competitive balance between usual radiative 
cooling processes and stochastic acceleration is a relativistic Maxwellian function also called 
'pile-up' distribution (Schlickeiser 1985; Aharonian, Atoyan, & Nahapetian 1986; Henri & 
Pelletier 1991), 



and where 7 max is simply the value of the individual Lorentz factor of the particles for which 
the acceleration time is equal to the cooling time. It corresponds to an energy distribution 
function of particles homogeneous and isotropic in the momentum space with a exponential 
cut-off at 7 m ax- Note that in the case of a power-law distribution function (with spectral 
index s > 2), the enthalpy of the plasma is dominated by lower bound of the particle energy 
range 7 m i n . For a 'pile-up', particles are mostly concentrated near 7 max and the dynamics of 
the plasma is mainly controlled by the high energy particles. 

The inclusion of an escape term will modify the above solution. The model presented here 
will break down if the escape time is much smaller or much larger than the characteristic 
acceleration time at the critical Lorentz factor 7 max . In the first case, acceleration will be 
much slower than the escape and no relativistic pile-up can be formed. In the second case, 
the relativistic particles will remain a long time before escaping (and cooling) and the emis- 
sion of the acceleration zone will be important. In the following, we exclude these two cases 
and we assume that the escape time is comparable to other times at 7 max , neglecting the 
emission of the acceleration zone. A proper inclusion of the escape term would modify the 
solution of the type given by equation (1), but the general shape would be the same ; a 
low energy part behaving like 7 2 when the acceleration/diffusion is very fast, followed by an 
energy cut-off. For sake of simplicity, we will thus use equation (1) and replace the 3 — r 
exponent by 1. The shape of the SED high energy tail is only weakly dependent on this 
approximation, and not strongly constrained by the observations. 




(1) 



2.2. The cooling zone 



In order to obtain the energy spectrum of emitting particles, we assume that the particles 
are accelerated as previously described in some localized region and are injected during some 
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time in a spherical zone where they cool freely. In this zone, we consider the standard kinetic 
equation in the continuous loss approximation with no escape term. It gives the evolution of 
the differential energy density of the particles ^±(7; t) with a Lorentz factor between 7 and 
7 + d7, 

fo n ±(r, t) + ^7(7; t)n ± ( r , t) = Q( r , t). (2) 

Particles source term Q{^;t) will include in fact both the fresh particles injection term 
Qmj(7;£) and the production rate Q pTO d(l',t) due to the pair creation via photon-photon 
annihilation, which will be developed in section 2.4. We take the following approximate 
form for the injection term : 

~ , ,x fn 7 2 exp(-7/7 max ) if ^ t< £ inj , 

<2mj(7;*) = < . (3) 

[0 otherwise, 

The factor 7(7; i) in the energy advective part of equation (3) is the continuous particle 
cooling rate. As mentioned above, charged particles can cool both via the synchrotron 
process or via the IC scattering of the previous synchrotron radiation field. We can thus 
write : 

7(7; = 7syn(7; t) + 710(7; *)■ ( 4 ) 



2.3. The radiative processes 

In the following we detail the equations used to compute the radiative processes. A 
'tilde' accent denotes a parameter expressed in the observer frame, otherwise in the blob 
frame. 



2.3.1. The synchrotron emission 

Assuming an isotropic particle distribution, the synchrotron cooling rate is given by the 
well-known formula: 

7syn(7^) = -hjnl 2 ] hyn = ~ U B , (5) 

O /7? gC 

where Ub = B 2 /8tt is the magnetic energy density. The synchrotron emission coefficient 
js(v) is obtained by performing the integration over the whole differential particle density 
of the mean emission coefficient for a single lepton averaged over an isotropic distribution of 
pitch angles Rcs( z ) (Crusius & Schlickeiser 1986; Ghisellini et al. 1988) 

jsfot) = ^| J d 7 n ± (7;t) J R cs W, (6) 
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Rcs(z) = 2z 2 



3z 



K A/3 (z)K 1/3 (z) - - {Kl /S (z) - K 2 /3 (z)) 



(7) 



with z = v/S^vb and = eB /2iim e c, K n being the McDonald function of order n. An 
accurate approximation of the function Rqs(z) is given in appendix A. 



2.3.2. The Inverse Compton emission 

In the same way, the Inverse Compton scattering cooling rate reads 

A fic(r,t) = J de 1 e 1 J de/C Joncs (ei,e,7)n syn (e;t), (8) 

where we consider the Compton kernel /Cj ones computed by Jones (Jones 1968; Blumenthal 
& Gould 1970) for an isotropic source of soft photons, considering the full Klein-Nishina 
cross section in the head-on approximation. More precisely, we have 

/C Jo nes(ei, 6, 7 ) = ^~f{q, T e )Q{q - 1/4 7 2 )6(1 - q), (9) 

Te = 46 7 , q = ^^ y (10) 

4e7(7 - ei) 



/(?,r e ) = 



2q]nq + (l + 2q)(l-q)+ 1 -^ q (l-q) 



(11) 



where G(i) is the usual Heaviside one-step-function and o"Th the Thomson cross section. We 
derive the IC emission coefficient Jq{v) in a similar way. The Compton kernel (equation 9) 
is this time integrated over the synchrotron emission spectrum : 

3c{yi\t) = ^e 1 J J /Cj oncs (ei,e, 7 )n syn (e;t)n ± ( 7 ;t)d 7 de, (12) 

where n syn (e) is the differential synchrotron photon density and €\ = hvxjm e (?. These 
equations are integrated numerically following the time evolution of the particle spectrum 
to find the time- dependent emission spectrum of the plasma. 



2.3.3. The pair creation process 



Gamma-rays photons produced in the blob can be absorbed by the photo- annihilation/pair 
creation process 7 + 7 — > e + + e~ (Gould & Schreder 1967a) for which the total cross section 
is 



a(x) 



3(7' 



Th , 



16 



1 -x 2 



(3 _ x 4 )ln i±^ _ 2x{ 2 _ x 2^ 

X Ob 



(13) 
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The 7-7 attenuation optical depth per unit length d£ reads 

d 



( tf r rr(ei) = J J J o-((3)n ph (e 2 ,{l)(l - ^)de 2 d 2 n, (14) 

where /3 = /3(ei,e 2 ,/i) = (1 — 2/e!e 2 (l — t 1 )) 1 ^ 2 is the velocity of the pairs in the center-of- 
mass system, €\ (resp. e 2 ) the energy of the high (resp. low) energy photon and \x the cosine 
of the collision angle (Coppi & Blandford 1990). For gamma-rays in the TeV range the pair 
production cross section is maximized when the soft photon energy is in the infrared range, 

EE h 

A = A ce 7 ~ 2.4 7 /im, where A ce = . (15) 

2m e c z 1 TeV m e c 

According to the previous relation, we can distinguish two different sources of soft photons 
being able to absorb the high energy tail of blazars. Firstly, through intrinsic attenuation, 
gamma-ray photons can interact with photons of the synchrotron bump in the source. As- 
suming again the synchrotron photon field to be isotropic in the blob frame, the integration 
over solid angle in equation (14) can be analytically computed. More precisely, one gets : 

d 1 f 

— r 77 (ei) = - / ^ P hx(e2)^ > pp(eie 2 ) de 2 , (16) 
where M pp (x) is the angle-averaged pair production rate (cm 3 /s) and reads 



■ pp {x) = c/ dfi - - - cr(P), Merit = max(-l, 1 - 2/x), 




B(x-l), 



:i7) 



and introducing the function 
1 



+ 



2u(2 + u) u 2u , . ,11 

/ + + 2 In 1 + u) - - + — 

(u + l) 2 4 u + 1 2 4u 



\n(u)+ 



7/12 7T 2 

+ 2dilog( M + l)-^ + ---^ + l + - (18) 
2 2u 1 + u 6 

and the dilogarithm function as (Abramowitz & Stegun 1964; Gould & Schreder 1967a) 



/X 
du 



(«-l) 



-^--iln 2 (x-l) + X;^-(x-l)- fc , x>l. 

k=l 
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To compute the photon escape probability ^ e sc(^i; t) of a photon of energy e\ = hvi/m e c 2 > 
1 (or also the the spectrum attenuation coefficient ^.bs^i! t)), we use the following approx- 
imate expression (Marcowith, Henri & Pelletier 1995), 

(19) 

with r 77 ~ i?dr 77 /d£. The factor between brackets is the usual solution of the transfer 
equation in the plane-parallel geometry approximation, and can approximate the photon 
escape probability in a spherical source of size R. The extra exponential factor in equation 
(19) has been introduced by Marcowith, Henri & Pelletier (1995) to account for the possibility 
for high energy photons to annihilate outside the source, because the soft target photons are 
not confined in the source like the particles; rather their density decreases slowly on typical 
length scale equal to the source radius. 

Secondly, the high energy photons can also interact with the photon field of the Diffuse 
Infrared Background (DIrB) (also called CIB for Cosmic Infrared Background) during their 
travel through the universe from the source to the observer (Gould & Schreder 1967b). 
Hereafter, we call this effect extrinsic absorption. DIrB is the extra-galactic light from the 
optical to sub-millimeter range, which records basically the history of star formation (for 
a review see Hauser & Dwek 2001, and references therein). If we ignore the secondary 
gamma-ray emission in the direction of the observer, it results that the emitted differential 
flux is attenuated by a factor (Gould & Schreder 1967b; Stecker, de Jager & Salamon 1992; 
Vassiliev 2000) 

Kb* (v, z.) = exp [-r;f ( Cl , z s )} , e, = hu/m c c 2 . (20) 
For close sources (z s « 1), the expression (14) gives: 

r^(e 1 ,z s ) ~ ^ f dp ^ [°° de 2 a(P)n DIrB (e 2 ), (21) 
J-i 2 J eth 

where e t h = 2/(1 — /x)ei, H is the Hubble parameter (assumed equal to 65 kms -1 Mpc -1 
throughout this paper), and n D i rB the density of the DIrB photon field. We have estimated 
this density by performing a Chebyshev interpolation using the measurements data compiled 
by Hauser & Dwek (2001) excluding the two points of COBE-DIRBE at 60 /im and 100 /tm. 
But note that with or without these points, the resulting absorption coefficients are quite 
similar above 9 TeV as shown in figure 1 and do not change the main results of this work. 

2.4. Pair production rate and pair cascade 

As mentioned above, the kinetic equation source term includes the contribution of the 
population of created particles in the pair production process as calculated above. Remarking 



1 — e T " n 



'77 



- 9- 



that a hard photon with a reduced energy e > 1 interacts preferentially with a soft photon 
of energy ~ 1/e < 1 to form a pair e + /e~ close to the production threshold. Consecutively, 
both particles have thus a similar energy 7 and we formally write the conservation of energy 
as e+ 1/e « e = 27. Then the pair production rate reads 



<2prod(7; t) 



dN 



dtd-fdV 



2n a bs0 = 2 7;^)- 



(22) 



Assuming the IC emission is isotropic in the plasma rest frame, the differential photon 
absorption rate density per energy and time unit, is given by 



n a bs(e) = 47rm e c 2 ^j^ ab s(^; t), 



(23) 



and we finally obtain 



Q P iod(r,t) 



iirm P c 



h 2 



^abs(^t) 



u=2'ym e c 2 / h 



(24) 



2.5. Time-averaged spectra 

At time t, the whole specific intensity in the plasma rest frame reads, 

S„(v; t)~R \j s {y- 1) + 3c (u; t) • <rf abs (z/; t)} ■ <g£{y, z s ), (25) 

where all parameters and emission coefficients are expressed above. All the physical quanti- 
ties must be converted from the blob frame to the observer frame, taking into account the 
Doppler boosting effect and the cosmological corrections according 

#*(£; 1) ~ 7r^U|(l + z a )S v {y\ t) with v = z>, (26) 

where di is the usual luminosity distance, 5b = l/r^(l — /3gcos0) is the Doppler beaming 
factor of the source and 9 the viewing angle. 

The observed spectrum is finally obtained by assuming that an observation takes place in 
the interval [t bs, ^obs + At Q b s ] (time t = is related to the beginning of the injection of fresh 
particles). The time-averaged spectrum is then 

n hs W = (Mv))i — ~~~ / M*\ t) at. (27) 
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3. Results and discussion 

3.1. General behavior 

Our model requires eight parameters. Three of them are related to the properties of 
the source, namely the magnetic field strength B , the radius R (or i? 15 when expressed 
in unit of 10 15 cm) and the bulk Doppler factor 5b- Three others characterize the injected 
plasma: they are the characteristic Lorentz factor of the 'pile-up' EDF 7 max , the number 
of injected particles which can be characterized by the integrated Thomson optical depth 
t xh = R&Th J dt J d7Q in j(7; t) and the injection time t in j. The remaining parameters are the 
observational ones, t Q bs and At G b s - 

In a steady-state model, there would only be five parameters since the last three ones would 
be irrelevant. More exactly the integrated Thomson optical depth r Th should be replaced 
by the constant particle optical depth r^ h = nRaTh, where n is the particle density in 
the source. For a 'pile-up' distribution, and neglecting pair creation, the whole spectrum 
is entirely characterized by two peak energies and two corresponding fluxes, corresponding 
respectively to the synchrotron and the IC bumps. Thus there would be only one free 
parameter left, which can be taken for instance as the unknown Doppler factor 5b ■ Exactly 
the same spectrum would be obtained by varying 5b and adjusting the other parameters 
accordingly. Since the TeV emission is dominated by the Klein-Nishina cut-off, the IC peak 
energy is simply proportional to 5s7 max , whereas the synchrotron peak energy is proportional 
to 5s7max-^- Thus the following scaling laws would apply : 

7max OC 5 B , , . 

B oc 5 B . (28) 

The integrated synchrotron luminosity scales as N^ &yi B 2 5 B) where N = AnR 3 n/3 = Attt^R 2 / (3o"Th) 
is the total number of particles in the source. So one gets the other following scaling law : 

R 2 r^ h oc 5 B \ (29) 

The final condition must be determined by the fact that the IC luminosity L IC is directly 
related to the synchrotron photon density in the source, which is itself directly related to 
the radius of the source (the synchrotron luminosity does not depend on this radius for a 
given number of particles). The magnetic energy density scales as B 2 so as 5 2 B (equation 
(28)), whereas the synchrotron photon energy density scales as Li C 5 B 4 R~ 2 . If one neglects 
the Klein-Nishina correction, the ratio of synchrotron to Compton luminosity is simply the 
ratio of magnetic to soft photon densityand is fixed by the observations. One would thus 
expect the following scaling law 

R oc 5 B 3 . (30) 
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In fact the real condition is more involved because the Klein-Nishina cut-off diminishes 
strongly the number of photons effectively available for IC scattering, in a way depending on 
the Doppler factor and the shape of the synchrotron spectrum. But qualitatively, one can 
always choose the radius of the source to adjust the IC luminosity to the observed value. 
There is however a limitation to the Doppler factor due to the existence of pair creation 
process. For to small values of 5b, the 7 — 7 optical depth, can increase so much that it 
becomes of the order of unity. Then the IC luminosity stops increasing, and rather it starts 
to decrease with decreasing radius because of 7 — 7 absorption. For a given Doppler factor, 
there is thus a maximum reachable IC luminosity. Conversely for a given IC luminosity, there 
is a minimum Doppler factor (which can be of course 1 in some cases where pair creation 
is never important). Note that the variability time scale is also a possible limitation, rapid 
variability requiring also high Lorentz factors. 

In principle the time- dependent model is more constrained. It requires three more free 
parameters (the injection time and the two observation times), but the entire shape of the 
synchrotron spectrum is depending on these times. One can see that by varying 5b, the 
cooling time of particles with the typical energy 7 max varies as (7max-B 2 ) -1 oc 5^ in the 
blob frame, so as 5g 2 in the observer frame. The shape of the synchrotron spectrum will 
remain unchanged if one scales all times proportionally to S^ 1 in the blob frame, or S^ 2 in 
the observer frame. However, one of these times, namely the observation lasting time At D b s 
is not a free parameter. Thus if the model fits perfectly well the data for all values of 5b, only 
one of these values is compatible with the actual value of At obs . So theoretically a unique 
set of parameters (if any) can fit the data. Of course things are not so ideal : because of 
observation error bars, data will be fitted by a set of possible values with a satisfactory x 2 
test. 



3.2. Approximated analytical solution of the kinetic equation 

For illustrative purposes, we develop here the simplest case where one can neglect the 
Inverse Compton cooling in comparison with the synchrotron cooling, and where pair pro- 
duction is unimportant. We can express analytically the general solution of equation (2) 
which satisfies the boundary condition n±( , y;t = 0) = for any arbitrary injection term 



Qinj(r,t) by, 




(31) 
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where t(7o, 7) is the energy drift time, i.e. the time needed for a particle of energy 70 to 
cool down to the energy 7, 

«™ ) = £wy (32) 

For the injection term chosen as in equation (3), equation (31) can be analytically integrated 
and we finally obtain, 



n±(r,t) 



•'•'o7max T-p f 3 a (r,t) \ _ T f 3 Kr,t) 



^syn7 2 



7max / \ 7) 



(33) 



where T(-,x) denotes to the incomplete gamma function. 

Here the parameters 0(7; t) and 6(7; t) represent respectively the lower and upper bounds of 
the integral (31), where the integrand does not vanish. To evaluate them, we distinguish three 
time intervals for each value of 7. Let us define the parameter t coo i(7) = | d'y / d'y | = l/yfcsyn- 
Note that t CO oi(y) = ^(00,7), i.e. it represents also the time spent by an initial infinite 
energy particle to cool down to 7. 

1. Initial stage where t < min(t cool (7), t in j) 

h{l]t) = T^-^t ; «r,t) = r 

Particles are still injected at the energy 7 but high energy particles have not yet cooled down 
to 7. Particles injected between 7 and some finite upper bound contribute to the integral. 

2. Cooling stage where max(£ inj , £0001(7)) < t < £ in j + £0001(7) 

7 

6(7; -»■ 00 ; a (r, t) = 



1 ksynlit ^inj) 

Particles are no more injected but some high energy particles are still cooling down to 7. 
Particles injected above a finite energy (larger than 7) contribute to the integral. 
3. Intermediate stage. For intermediate values of t, we must distinguish two specific energy 
ranges. We define a critical value of the individual Lorentz factor of the particles, 7i im = 

1 / k S ynt i n j . 

3. a) low energy range where 7 < 7 lim (or t inj < t cool (y)), 

Kr,t) = -, — ? — i ; a (r,t) 



1 ^syn7^ 1 ^syn7(^ ^inj) 

Injection is finished but very high energy particles have not yet cooled down to 7. Particles 
injected in some interval above 7 contribute to the integral. 

3.b) high energy range where 7 ^ 7u m (or t in j ^ £0001(7)), 



-13- 



6(75*) ^oo ; 0(7; t)=7- 

Conversely, very high energy particles have time to cool down to 7 while injection of fresh 
particles still takes place. All particles injected above 7 contribute. Note that this is the 
only stage for which 72(7; t) does not depend on time. A steady-state is set during this stage 
(although not in the spectrum because the whole distribution is not steady). 

We could also define an end stage for which t > £0001(7) + tmj and where ^±(7) = 0. 

Introducing the reduced variables, e = 7/7max and = tj/t C ooi(7max), the previous 
equations can be collected into the following relation, 

n( £ ;r) = ^%^ inj ( £ ;r), (34) 

where 



uj Tin . (e; r) = e 2 0(1 - e max(0, r - r inj )) x 



x 



r 3, 



1 — e max(0, r — r in j) 



0(1 -er)r 3, 



1-er 



• (35) 



An example of a resulting cooling pair-EDFs at different times is plotted in figure 2. One 
clearly sees the initial stage where the EDF is built, the formation of a 7 -2 EDF due to the 
cooling and the subsequent cooling of the whole distribution after the injection has stopped. 
As we will see in realistic simulations, the shape is however strongly modified when taking 
into account the IC cooling process (which is not simply dependent on the energy because 
of the Klein-Nishina cut-off) and the pair production term. 



3.3. Application to Mrk 501 data 

3.3.1. Observations 

We have applied the model to fit the spectral energy distribution of Mrk 501 during the 
period of 1997 April when this source experienced an intense period of activity. From this 
period, we distinguish two different activity states, namely the 'high state' from the April 16 
and the 'medium state' form the April 7. Simultaneous data are taken from BeppoSAX for 
the X-ray observation (Pian et al. 1998) and from the French Atmospheric Cerenkov Tele- 
scope CAT for the spectrum in the TeV energy regime (Djannati-Atai et al. 1999; Barrau 
et al. 1998). In a first step, we have corrected the high energy spectra using the attenuation 
coefficient computed previously. Note that the last corrected data point of the high state 
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may be not meaningful, leading to a concave up corrected high energy spectrum. The most 
simply and obvious explanation of this problem is an over estimation of the measured high 
energy tail of the blazar spectrum or/and of the DIrB density itself. 

The value of the At G b s is constrained by the observing duration. However there are some 
discontinuities in the observing time, which make the global observation period length dif- 
ferent from the real exposure time. We have assumed that on average our At Q b s corresponds 
to the global observing time of the BeppoSAX instruments LECS/MECS and PDS (~40000 
s for the April 16, 37845 s for the April 7). 



3.3.2. Solutions without pair creation 

We can reproduce the data with parameters corresponding to almost no pair creation. 
The parameters used to fit the data are reported in Table 1 and the resulting synthetic SED 
are plotted in Figure 3. The fits are quite satisfactory for the high energy part of the SED. 
One-zone models are appropriate only to to reproduce this high energy, variable part. They 
can not account for the radio emission produced at much larger scale, where the whole jet 
contributes, possibly being the superposition of many successive flaring events. 
In some range of Lorentz factor, a steady-state solution corresponding to 7~ 2 power-law 
(see equation (33), case 3.b)) can be observed giving a z/ -1 / 2 synchrotron flux index. This 
corresponds to the main part of the April 16 spectrum and and the low energy part of 
the April 7 one. Above this range, the spectrum is modified by a factor £ coo iAobs because 
the particles have cooled before the end of the observation. This produces a steepening 
of the spectrum by Act = 0.5 and explains the flat high energy part of the synchrotron 
spectrum. The position of the spectral break is thus directly determined by the ratio between 
the cooling and the observation time. We obtain values of the magnetic field and of the 
transverse radius of the source (R15 ~ 1) in agreement with other models (see e.g. Bednarek 
& Protheroe 1999; Tavecchio et al. 2001; Ghisellini, Celotti, & Costamante 2002; Katarzyhski 
et al. 2001; Konopelko et al. 2003, ). It turns out that this minimum Doppler factor implied 
by 7 — 7 absorption in steady state models is quite high for Mrk 501, as noticed already 
by several authors. When the IC luminosity is corrected from extrinsic absorption, it can 
be as high as 50 for stationary models (see e.g. Bednarek & Protheroe 1999; Konopelko et 
al. 2003). In the time dependent model, the constraint is somewhat released because high 
energy photons (TeV) are not produced at the same time as low energy, IR photons that 
are mainly responsible for their absorption. So the actual density of soft photons during 
the emission of highest energy TeV photons is lower than that measured on average. We 
see that good fits can be obtained with Doppler factors around 25, which implies a bulk 
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Lorentz factor larger than 12. Not that this still relatively high value is necessary due to the 
correction of the spectra from from the DIrB absorption. These high values are generally in 
conflict with other features of TeV BL Lacs ; first VLBI/VLBA sub-parsec observations of 
Mrk 501 and others TeV blazars which don't exhibit superluminal velocities at mas scale or 
excess in the derived synchrotron brightness temperature (Edwards & Piner 2002; Piner & 
Edwards 2004). Moreover, if we suppose that BL Lacs are the beamed counterparts of the 
FR-I radio galaxies, the ratio of their bolometric luminosity should be of the order of 8 B . 
A 8b larger than 10 would result in a much larger contrast than what is observed (Urry & 
Padovani 1995). The spatial density of beamed vs unbeamed objects seems also to be better 
consistent with a Lorentz factor of around 3. All these argument disfavour a high Lorentz 
factor. This is in fact a problem inherent to any one zone model. We will argue in a future 
work (Sauge & Henri in preparation) than inhomogeneous models may solve this issue. 

Like B and 8b, 7max is kept constant between the two different states. In the light of 
our scenario, the 'medium state' spectrum could be due to a previous ejection observed in a 
later stage (with respect of the injection time) than the April 16 one and in a much larger 
part of the jet. 

We can also estimate the minimum variability time scale deriving from standard causal- 
ity arguments and given by t var > t V ar,min = 555 -R15 (1 + z s )8 B ' 1 min. Its value is reported in 
the last column of the Table 1 and is equal to 15 min for the high state and roughly 40 min 
for the medium state where the size of the source is much larger as mentioned above. Note 
however that the injection time is much larger than the above values. 

Simulated light curves corresponding to April 16 parameters are shown in figure 4. 
The curves are calculated for the 3 energy ranges of BeppoSAX and for very high energy 
(VHE) instruments above 250 GeV. Lags between various energy bands are clearly visible 
on the curves. Remarkably, the high energy photons lead the soft energy ones in the X- 
ray synchrotron component, but the high energy gamma-ray curve lags the synchrotron 
component. This is due to the fact that when particle distribution cools, the photon density 
below the Klein-Nishina limit increases. This makes the IC emission level keep rising even 
after the injection has stopped. Precise comparisons with observed light curves have not 
been made because the statistics is too poor to allow a meaningful analysis. However we can 
see qualitatively that complex temporal effects can arise from a time-dependent simulation 
leading to apparently contradictory behaviours. Note that whereas the presented simulations 
span the whole temporal interval of BeppoSAX observations, high energy observations have 
been performed only during a limited time within this interval. If they have taken place 
around the maximum of the light curve, few variations are expected. 

Of course, a more realistic model could include a succession of different flares as well as 
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a variation of the injection characteristic energy. This is however beyond the scope of this 
paper which only aims at showing that 'pile-up' particle energy distribution can well explain 
the high energy shape of TeV blazar spectra. 

3.3.3. Solutions with pair creation 

Interestingly, we found also possible solutions with a much lower bulk Doppler factor, 
and accordingly a lower magnetic field, a higher particle Lorentz factor and a larger radius. 
An example of such parameters is shown in Table 2, and corresponding spectrum is displayed 
in figure 5. In this regime, pair production can be important. Although the IC flux should 
be much larger, the 7 — 7 absorption reduces effectively the luminosity so that it can be 
compatible with the observations. The intense creation of new particles produces more 
synchrotron photons, which accelerates the cooling. This phenomenon also amplifies the 
effect of the increase of the Klein-Nishina threshold energy. This leads eventually to a 
catastrophic pair production/cooling process, producing a strong flare at GeV energy and 
in low energy X-rays. It explains the bump occuring in the GeV/sub-TeV energy range 
(and to a lesser extent in the radio sub-millimeter range). The flare is also more clearly 
visible on the light curves (figure 6) appearing as a very sharp flare in some energy ranges. 
Although the relevance of such scenario is not clear, it may be possible that such events 
could explain the most rapid variations in the TeV light curves. The GeV flare would have 
been in principle observable by the EGRET instrument. However due to the small number of 
simultaneous EGRET/TeV observations and the briefness of the event, it could have escaped 
any detection. We note also that for different value of the parameters, it is possible that 
the X-ray flare shifts toward lower energy, disappearing from the X-ray data. This could be 
related to the 'orphan flare events' observed in some occasions (Krawczynski et al. 2004). 

4. Conclusion 

This paper shows that the high energy spectra of TeV blazars can be well reproduced by 
a cooling 'pile-up' EDF. This offers an alternative to the narrow power-law injection terms 
often used in the literature, whose justification is unclear. Inclusion of time-dependent effects 
permits to reproduced the main features of both the light curves and the time-averaged 
spectra. For a given SED shape, the parameters of the model are fully constrained. For 
Mrk 501, different states could be the result of the time variation of the transverse radius 
of the source R, the quantity of injected leptons via r^h and the observational parameters 
with respect of the injection time. However, this one zone model shares a common issue 
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with other homogeneous models : it requires high values of the bulk Lorentz factor to avoid 
a strong gamma-ray absorption, even in the case of the optically thick solution (see section 
3.3.3). These high Lorentz factors appear to be inconsistent with those deduced from FR 
I/BL Lacs unification models (Urry & Padovani 1995; Chiaberge et al. 2000). They are also 
difficult to reconcile with the absence of superluminal motion at TeV scale and relatively low 
brightness temprature of TeV blazars (Edwards & Piner 2002; Piner & Edwards 2004). We 
will argue in a future work that inhomogeneous models could fix this issue. 

We thank A. Djannati-Atai, P. Espigat, F. Piron, S. Pita and A. Lemiere from the CAT 
collaboration and A. Marcowith for helpful interaction. Critical remarks of an anonymous 
referee helped to improve significantly the final version of this paper. L.S. would like to 
thank J. Ferreira and P.-O. Petrucci for very enthusiastic and stimulating discussions and 
encouragement. Part of the simulations reported here has been performed at the " Centre de 
Calcul Intensif de VObservatoire de Grenoble.'' 1 . 



A. An accurate analytic approximation of pitch angle averaged 

synchrotron emitted power 

As mentioned in section 2.3.1, the Rq^{x) function (see equation7) results from the 
monochromatic emitted power averaged on a population of particles with randomly dis- 
tributed pitch angle. It mathematically reads (Crusius & Schlickeiser 1986) 

Rcs(x) = \ r d9 sin 2 6F syn (-?-) , (Al) 
z Jo vsmc?/ 

where F syn (x) is the usual synchrotron fundamental kernel (Blumenthal & Gould 1970), 

POD 

F syn (x)=x / dzK 5/3 (z). (A2) 

J X 

For i < 1 approximate expression of F syn is 

F syn (x)n J* — (^y /3 , (A3) 



V3r(l/3) V2 

and we immediately obtain the relevant one for -Res, 

2 l/3 



R^( x ) = - d9 sin 2 6F syn (-—) = —T 2 {l/3)x 1/3 w 1.808418 x 1/3 . (A4) 
z Jq v sin u/ o 
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Conversely, for large argument (x >> 1), the asymptotic development of F syn is (nx/2) ' e x 
and then 

= \J\e *?9F„ (-^) = \ (f ) 1 7*d ( - sin V»« (A5) 

The preceding integral can re-writes as, 

l(x)= T d9exp(-f x (9)) , = ^-^Ln(sin0). (A6) 

Remarking that the exponential argument is symmetric around 6 = it/ 2 and therefore on 
the integration range, we set 6 = n/2 + and make a taylor expansion around ip 

= ^Tn ~ \ ln ( C ° S ™ X +'t( X + D + (A7) 

COS if 2 2 \ 2 J 

Then integral (A6) rewrites, 




which is a standard gaussian integral. Noting that larger is x sharper is the integral argument, 
lower and upper bound can be extended to infinity as x ^> 1 (error is less than 10~ 10 for 
x ^ 3 and ). We finally obtain, 

*» w - f ^fr~'- (A9) 

For intermediate values, we perform a polynomial fit of the form, 

flint 0*0 = a (x ai + a 2 x 2ai + a 3 x 3ai ) exp(-a 4 x a5 ), (A10) 

where the coefficient is obtained from least square fitting and is given in table 2. In this 
regime, RMS error is in order of 0.05% . 

REFERENCES 

Abramowitz, M., & Stegun, LA. 1964, Handbook of Mathematical Functions, Applied Math- 
ematics Series - 55 (1965, New York: Dover Publication Inc.) 

Aharonian, F. A., Atoyan, A. M., & Nahapetian, A. 1986, A&A, 162, LI 

Aharonian, F. et al. 2002, A&A, 384, L23 



-19- 

Aharonian, F. et al. 2003, A&A, 406, L9 

Barrau A. et al. 1998, Nucl. Instrum. Methods Phys. Res A, 416, 278 

Bednarek, W. & Protheroe, R. J. 1999, MNRAS, 310, 577 

Blumenthal G., Gould R. 1970, Rev. Mod. Phy., vol 2, 2, 237 

Catanese, M., et al. 1998, ApJ, 501, 616 

Catanese, M. & Sambruna, R. M. 2000, ApJ, 534, L39 

Chadwick, P., et al. 1999, ApJ, 513, 161 

Chiaberge, M., et al. 2000, A&A, 358, 104 

Coppi, PS. & Blandford, R.D. 1990, MNRAS, 245, 453 

Crusius A., Schlickeiser R. 1986, A&A, 164, L6 

Djannati-Ata'i, A. et al. 1999, A&A, 350, 17 

Djannati-Atai, A. et al. 2002, A&A, 391, L25 

Dermer, C. D., Miller, J. A., & Li, H. 1996, ApJ, 456, 106 

Edwards, P. G. & Piner, B. G. 2002, ApJ, 579, L67 

Finkbeiner, D. P., Davis, M., & Schlegel, D. J. 2000, ApJ, 544, 81 

Georganopoulos, M. & Kazanas, D. 2003, ApJ, 594, L27 

Ghisellini, G., Maraschi, L., & Treves, A. 1985, A&A, 146, 204 

Ghisellini, G., Guilbert P., Svensson R. 1988, ApJ, 334, L5 

Ghisellini, G., Celotti, A., & Costamante, L. 2002, A&A, 386, 833 

Gould, R. J. & Schreder, G. P. 1967a, Physical Review , 155, 1404 

Gould, R. J. & Schreder, G. P. 1967b, Physical Review , 155, 1408 

Hauser M.G., Dwek, E. 2001, ARA&A, 39, 249 

Henri P., Pelletier G. 1991, ApJ, 383, L7 

Holder, J. et al. 2003, ApJ, 583, L9 



-20- 



Horan, D., et al. 2002, ApJ, 571, 753 
Jones F.C. 1968, Phys. Rev., 167, 1159 

Jones, T. W., O'dell, S. L., & Stein, W. A. 1974, ApJ, 192, 261 
Jones, F. C. 1994, ApJS, 90, 561 
Katarzynski K., Sol H., Kus A. 2001, A&A, 367, 809 
Konigl, A. 1981, ApJ, 243, 700 

Konopelko, A. K., Mastichiadis, A., Kirk, J. G., de Jager, O. C, & Stecker, F. W. 2003, 
preprint (astro-ph/0302049) 

Krawczynski, H., et al. 2004, ApJ, 601, 151, preprint (astro-ph/0310158) 

Lacombe, C. 1977, A&A, 54, 

Marcowith, A., Henri, G., & Pelletier, G. 1995, MNRAS, 277, 681 

Nishiyama, T., et al. 2000, AIP Conf. Proc. 516, Proceedings of the 26th International 
Cosmic Ray Conference, ed. B.L. Dingus, D.B. Kieda & M.H. Salamon (Salt Lake 
City, Utah: AIP), 3, 370 

Pelletier G. 1985, in proceedings "Plasma turbulence and Astrophysical Objects", SFP 

Pelletier G., Sol H. 1992, MNRAS, 254, 635 

Pian, E. et al. 1998, ApJ, 492, L17 

Piner, B. G. & Edwards, P. G. 2004, ApJ, 600, 115 

Punch, M., et al. 1992, Nature, 358, 477 

Protheroe, R. J. & Stanev, T. 1999, Astroparticle Physics, 10, 185 
Quinn, J., et al. 1996, ApJ, 456, L83 

Rybicki, G. B. & Lightman, A. P. 1979, Radiative Processes in Astrophysics (New York: 
Wiley-Interscience) 

Schlickeiser, R. 1985,A&A, 143, 431 

Stecker, F. W., de Jager, O. C, & Salamon, M. H. 1992, ApJ, 390, L49 



- 21 - 



Tavecchio, F. et al. 2001, ApJ, 554, 725 
Trussoni, E. el al. 2003, A&A, 403, 889 
Urry, C. M. & Padovani, P. 1995, PASP, 107, 803 
Vassiliev, V. V. 2000, Astroparticle Physics, 12, 217 



This preprint was prepared with the AAS IATgX macros v5.2. 



-22 - 



Table 1. Mkn 501 physical fit parameters 





B 


Rib 


8b 


7max 


T"Th 


^inj 




^obs 




^var,min 




(G) 


(10 15 cm) 




10 6 


io- 10 


(ks) 




(ks) 


(ks) 


(min) 


high state 


0.077 


0.65 


25 


1.26 


58.5 


9.4 


3.29 (: 


= 0.35 t inj ) 


39.7 


15 


medium state 


0.075 


1.75 


25 


1.29 


8.95 


40 


35 (= 


0.875 t inj ) 


37.8 


40 


'GeV flaring state' 


0.047 


1.06 


16 


2 


144.1 


24.4 


8.54 (: 


= 0.35 t inj ) 


102 


38 
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Table 2: Resulting fitting coefficients from the polynomial approximation of the function Rqs 



a = 0.201 447x10 

ai = 0.344 606 
a 2 = -0.429 682 

a 3 = 0.273 331 xl0~ 2 

a 4 = 0.966 844 

a 5 = 0.964 518 

X 2 = 1.27xl0~ 5 

RMS % error = 0.05 
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Fig. 1. — Synthetic absorption coefficients calculated for Mrk 501 (z s ~ 0.034, Hq = 
65 kms" 1 Mpc -1 ). Solid curve (high hypothesis) takes into account the 100 and 60 /im 
flux reported where dashed curve (low hypothesis) does not. 
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Log 10 [ E=y/ Ymax ] 

Fig. 2. — An example of resulting cooling 'pile-up' energy distribution function w Tinj (e; r) 
for Tinj = 4 a time r = 1,2,3,4,5. Also represented, in dashed straight line, a power-law of 
index 2 typical which results from the radiative cooling of a mono-energetic particle energy 
distribution. 
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Fig. 3. — Fits of Mrk 501 data during the flaring period in 1997 April by the model. See 
Table 1 for the detail of the parameters values. For each state, the solid curves show the 
emergent spectrum modified by the intrinsic absorption due to the pair creation process 
in the spherical blob. The dashed curves show the unabsorbed one. Filled gray circles : 
observed CAT data. Open squares : data corrected from the DIrB attenuation. 




Fig. 4. — Mrk 501 April 16 simulated light curves for BeppoSAX instruments (left panel) 
and above 250 GeV (right pannel). Upper axis is graduated in unit of hours in the observer 
frame. 
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Fig. 5. — Time-averaged homogeneous one- zone SSC model for spectra of Mrk 501 in the case 
of an intense pair production. Thick grey lines show the EGRET and GLAST instruments 
sensitivity curve (3a during 102 ks). 
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Fig. 6 



- Same as figure 4 for the case of an intense pair production. 



